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ABSTRACT 

Delaunay- Voronoi mesh systems provide a generalization of the classical rectangular stag- 
gered meshes to unstructured meshes. In this work, it is shown how such “covolume” dis- 
cretizations may be applied to div-curl systems in three dimensions. Error estimates are 
proved and confirmed by a numerical illustration. 
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1 Introduction. This paper deals with the problem of numerically determining 
a three dimensional vector field from its divergence and curl. The need to do this 
occurs in several fields including fluid dynamics and electromagnetics. The governing 
equations are often solved using indirect approaches including potential formulations 
or Biot-Savart type integrals [5] and least squares [1]. 

A good reason to be cautious about direct discretizations is that the equations are 
overdetermined, being a system of four equations in three unknowns, and it is not 
immediately evident what effect this will have at the discrete level. Nevertheless, 
since potential formulations can suffer from spurious mode problems [6] and the 
Biot-Savart approach requires special handling of boundaries, a direct treatment 
may be preferable. We present such a treatment in this paper. 

The method we discuss is a generalization of that presented in [3] for planar systems. 
As in the planar case we use dual mesh systems made up of “complementary vol- 
umes'' (or “covolumes” ) of which the prime examples are classical Voronoi-Delaunay 
mesh pairs. In the Voronoi-Delaunay approach the domain is suitably partitioned 
into tetrahedra which are considered to form a primal (Delaunay) mesh in the do- 
main Q. A dual (Voronoi) mesh is formed by connecting the circumcenters of the 
tetrahedra, giving a set of convex polygons each of which encloses a unique tetrahe- 
dral vertex. Combinatorial and topological properties of the primal and dual mesh 
systems play an important part in our analysis, allowing us to conclude, for example, 
that the solution of the discretized equations is exactly determined by the data if the 
same is true of the original system. On the other hand, the metric properties of the 
mesh system are used in deriving the error estimates. The interplay of topological 
and metric properties is a common feature of covolume discretizations. 

The paper contains six sections. The second section deals with combinatorial prop- 
erties of the mesh system. The third section shows how discrete vector fields are 
defined on the meshes and defines the basic integral operations for the discrete 
fields. These definitions are set up so that an important orthogonality property of 
vector fields is preserved (theorem 4). The discrete approximations and the basic 
error estimates are in section 4. In section 5 we specialize the results to rectangu- 
lar meshes for which we are able to obtain improved convergence rates. The final 
section contains some basic numerical verifications. 


2 Mesh notations and some basic results. We will prove our results for the 
case of a domain Q which is a bounded polyhedron in R 3 with boundary T. The 
results we present are generally valid for multiply connected domains and domains 
with cutouts. The techniques for extending the results are similar to those used in [3] 
for the tw r o dimensional case although there is greater variety in three dimensional 
cases depending on the genus of the surface(s) bounding the obstacle(s). 
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Assume that fi has a primal family of finite element style tetrahedral partitions, 
parameterized by the maximum side length which is generically denoted by h. We 
wiD assume that the ratios of the radii of the circumscribing spheres and the in- 
scribed spheres of all the individual tetrahedra are uniformly bounded above and 
below as h approaches 0. A dual mesh is formed by connecting adjacent tetrahedral 
circumcenters and, in the case of tetrahedra with a face on a boundary, by connect- 
ing their circumcenters with those of their boundary faces. By elementary geometry 
the connecting edges are perpendicular to the associated tetrahedral faces. These 
connections also form the edges of a set of polyhedra. It follows from elementary 
geometry that the edges of the tetrahedra are perpendicular to and in one-to-one 
correspondence with the faces of the dual polyhedra (covolumes) and dually. The 
reciprocal orthogonality between edges and faces is the key to the results which 
follow. A survey of algorithms for generating Voronoi-Delaunay meshes is given 
in [4]. 

For an entirely arbitrary tetrahedral partition the structure of the covolume partition 
can be complicated. To obtain a practical algorithm we will restrict the tetrahedral 
partition so that the dual partition has two (possibly dependent) properties. These 
properties are (1) interior covolumes have faces which are simple (planar) polygons 
and (2) each tetrahedral interior node lies in exactly one covolume. This does restrict 
the primal tesselation somewhat but not unduly. For instance the two conditions 
are certainly satisfied for a tetrahedral partition in which the dihedral angles are all 
acute. In that case, the two meshes form a Delaunay-Voronoi pair and the interior 
covolumes are convex. In addition the faces of the covolumes are not merely simple 
but convex as well. More generally any Delaunay-Voronoi pair will satisfy (1) and 
( 2 ). 

The N nodes of the tetrahedral mesh are assumed to be numbered sequentially 
in some convenient way, and likewise the T nodes of the dual mesh. Similarly, 
the F faces (edges) and E edges (faces) of the primal (dual) mesh are sequentially 
numbered. A subscript 1 on one of these quantities denotes the corresponding 
interior number. For example F\ denotes the number of tetrahedral interior faces. 
The individual tetrahedra, faces, edges and nodes of the primal mesh are denoted 
by Ti, pj, <Jk and i respectively. Those of the covolume mesh are denoted by primed 
quantities such as o'- . A direction is assigned to each primal edge by the rule that 
the positive direction is from low to high node number. The dual edges are directed 
by the corresponding rule. 

We will use the following standard results: 

Lemma 1 Let N , F, E and T denote the number of nodes , faces, edges and tetra- 
hedra of a given tesselation of a polyhedral domain. Then 

N + F = E + T + 1 
N\ + Fi = Ei + T - 1. 
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( 1 ) 

(2) 



These can easily be proven by successive deletions of tetrahedra from the triangu- 
lation. ■ 


3 Discrete vector fields. The main theorems in this section (theorems 3 and 
4) are discrete analogs of theorems in vector field theory. Theorem 3 asserts the 
existence of a scalar potential when discrete circulation equations are satisfied and 
theorem 4 gives a discrete Helmholtz decomposition of a field of normal components. 

Each boundary df.ij ( dfi ! k ) is assumed to be oriented by a right hand rule applied 
relative to the directed edge crj 

For each strictly interior dual edge a j we can form a vector whose kth component is 
the sign of the orientation of the edge relative to the orientation of the kth strictly 
interior dual face. From these vectors we obtain the F\ x E\ matrix G defined as 
follows: 

{ 4*1 if CTj is oriented positively along dfi' k 
— 1 if (Tj is oriented negatively along d\i\ 

0 if a j does not meet dfi' k . 


G has rank E\ - N\. To see the reason for this, note first that each column of G is 
associated with a dual face. Now take any interior covolume and sum the associated 
columns, each weighted ±1 according to whether its normal points out of or into the 
covolume. Then the contributions from each dual edge appear twice in each sum 
with opposite signs. In this way we can form one null vector of G l for each interior 
covolume. It is clear that there are no other vectors in the null space so the result 
follows. 


In addition to G we will use another matrix B\ containing the orientations of the 
tetrahedral faces relative to the outer normals of the tetrahedra. This matrix is 
defined for interior faces and has dimensions F\ x T. The definition of B\ is: 


{Bxh 


+ 1 if fij is oriented with dr x 
— 1 if f i 3 is oriented against dr x , 
0 if fXj does not meet dri . 


where “oriented with” means that the normal direction of fij is parallel to the outer 
normal of dr{. The T dimensional vector with 1 in every position is in the null space 
of B\ and is the only such vector, so that the rank of B is T — 1. 

In addition to B\ we will also use the matrix B which is obtained from B\ by aug- 
menting it with rows corresponding to the dual edges which are normal to boundary 
faces. Recall that these edges pass through the circumcenters of the boundary faces. 
The additional values for the domain of B are associated with these circumcenters. 

B\ and G are related through G t B\ = 0. To see this, take a standard basis vector 
and multiply it by B. Then values ±1 or 0 result on the dual edges the positive 
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(negative) sign being associated with a dual edge directed out of (into) the node 
which carries the unit value. Multiplication by G 4 corresponds to forming a signed 
sum of these values around the dual faces where the signs are precisely such as to 
effect a cancellation of the nonzero contributions. Transposing this result shows 
that R(G ) C N(B\) where N(-) and R(-) denote the null space and range of their 
arguments. 

Theorem 1 Let v € R Fi ■ Then 3 <j> € R T such that v - B\<f) if G l v = 0. 

Proof. Since Bi is Fi x T and of rank T - 1 and using lemma 1 we have 

dim N(B[) = F\~T+\ 

= E l -N 1 . 

Recalling that R(G) C N{B\) and since dim R{G) = E\ - N\ we have N ( B\ ) = 
R(G). Solvability of the equation 


B\(j> = v 

holds if 

( v, z) = 0 £ A ( ), 

where ( , ) denotes the standard Euclidean inner product. From above this is equiv- 
alent to 

(v,G4>) = 0 


and these are equivalent to the theorem. ■ 

The covolume algorithm for the div-curl equations provides approximations to the 
quantities u-n, where u denotes the exact solution vector at the center of a primal 
face and n denotes a normal to the face. The corresponding approximate quantities 
are denoted by u 3 where j indexes the faces. 

Any set of normal components defined on faces can be identified with R F . We will 
introduce an inner product into R F by 

(u,v) w := u j l 'j s j h 'j- 

where sj is the area of and /?' is the length of cr'. The resulting inner product 
space is denoted by U. and Uq denotes the space 

U 0 := {u eU;u | r = 0}. 
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The associated norm is denoted by ||*|| vv . Clearly, this is three times a discrete L 2 
norm. 

For each tetrahedron 77 discrete flux and divergence operators are defined on U by 

( Du) { := Y. u 3*i 

and 

(Du), := (Du),/ A,. 

By Sj we mean sj negatively signed if the corresponding velocity component is 
directed towards the inside of 77 and positively signed otherwise. 

For each interior covolume face n* k discrete circulation and curl operators are defined 
by _ 

(Cu)k := 

a j£ d »'k 

and 

(Cu) k (Cu) k /s' h . 

The tilde on /*' means that this quantity is to be taken with a negative sign if the 
dual edge is directed against the positive sense of description of d{i f k and a positive 
sign otherwise. 

Denote by S and H b the diagonal matrices 5 := diag(s J ) and H b diag(ft'). It 
can be checked directly that 

D = B l S 
C = G l H\ 

where H ; denotes the restriction of H f b to interior dual edges. We also define differ- 
ence operators 

Pb 4 > ■■= (H' b )~ x B<t> V<j> € R t 
P< t> := (H')~ x B\<f> V<peR T . 

Corresponding to theorem 1 we have 

Theorem 2 Let v € R Fl . Then 3 4 > € R T such that v = Pcf> if Cv = 0. ■ 

Now introduce two subspaces of Uq\ 

Z 0 = {tie U 0 \ Du = 0} 

Wo = {«£ Uq\Cu = 0}. 


Then we have 
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Theorem 3 Let u € Zq, v € Wo, then 

(u, v) w — 0. 

Proof. From theorem 2 and using summation by parts we have 

(u.v ) w = (u,P4>) w — {u,P b <f>) w = (Du,4>) V<£ € R T ■ 

Since u € 2 0 we have jDu = 0 and the theorem follows. ■ 


4 Div-curl systems. In this section we formulate the discrete approximations 
for div-curl equations and show that they have a unique solution (theorem 5). The 
basic error estimate is given in theorem 6. 

We consider the three dimensional div-curl problem for u, 



div u - p 

(3) 


curl u = u) 

(4) 


u n | r = / • 

(5) 

We will assume that the following compatibility conditions hold: 



o 

II 

3 

> 

(6) 


f pdv — f fds . 

Jn Jr 

(7) 

See [2] for the mathematical background to this problem. 
Equations (3)-(5) are discretized by 



II 

(8) 


II 

O 

(9) 


« lr = /> 

(10) 

where 

II 

JH ^ 

a. 


Pi 


(£y)* 

=r A* / U-tds 
S k 


f 

= — f fds pj e r. 

S 3 J»3 



The equations (8)-(10) are not all independent. In fact, there are at least E x relations 
between the equations (9) corresponding to the fact that the area weighted sums of 
the data over covolumes is zero (corresponding to the compatibility condition (7)) 
as well as a further relation which is equivalent to (7). The discrete equations do 
have a unique solution as we prove next. 
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Theorem 4 Equations (8)-(10) have a unique solution. 


Proof. Consider the homogeneous equations with = 0 and / = 0 in (8 )-( 10). 
These equations have at least one solution u, and this solution satisfies u G Z 0 and 
u G VV'o- But then theorem 3 implies u = 0 which implies the uniqueness part of the 
result. Let A denote the (T + F\ x F\ coefficient matrix of the linear system; then 
it follows that K has rank F \ . 

For the existence, we consider the linear system generated by (8)-(10). and let K t 
denote the corresponding augmented matrix of K of order (T + Ei) x ( Fi + 1). It 
follows from the compatibility conditions that this matrix has rank at most T - 1 + 
E\ — N\ = F\ and existence follows from this. ■ 

Next, we introduce two functions and defined by 

:= — / u*n ds fi 3 G f \ 

s j Jvj 

u[ 2 ^ •= tt f u n ds G 0. 

h j M 3 

If a'j is connected to the boundary T, := u^\ 


Lemma 2 Assume that u 6 Vr 1,p (TI),p > 2, then we have 

< c/l l u L P ,n- 


- u< 2 ) 


Proof. Let kj denote the polyhedron associated with the face and the dual edge 
<r' which is obtained by connecting the circumcenters of the primal cells which share 
fij to its vertices. By a Sobolev embedding theorem we have 

LUk)) p> 2 r = 1, 2, 

where /cj = cr' and k 2 = It follows that V : u -> - uf ] defines a bounded 

linear functional on W 1,p (Kj). Clearly, u' 1) - is zero for constant u so we have 

u( P | ^ Ci(fc,p)|u| liJVV 

and a standard rescaling argument then shows that the right side is 
where C is independent of h and u. Summing over the domain we obtain 

£ la'" - ufl v; < 



Then by Holder’s inequality and taking C such that h^F < C we have 

) 2/p / \ (p- 2 )/p 

E 1 

/ 

< Ch 5 - 6/p F( p - 2)/p |u£ p>n 

< C^ 2 l u li,p,n 

which we wanted to show. ■ 

From this lemma we can obtain the error basic estimate for the approximate solution: 


u' 1 * - u (2) 


= Ch 


5— 6/p 


Theorem 5 Assume that u € W l ' p (Q,),p> 2. ITien we have 


\U - U 


(Dl 


< C/i|u| lpn . 


Proof. Define r< ! > = u - u (l \i = 1,2, clearly, t’ (1) G 2o and u {2) G W 0 . Let 
u = (nf 1 * + u* 2 ))/ 2. A calculation using Theorem 3 shows that 


On the other hand, we have 

u - u - u - - (u^ - «( 1 ^)/2. 


Then 

and the result 


u - 


< P - «llw + 2 


111, ,(2) _ „(!) 


follows from lemma 2. 


w 


5 Uniform meshes. In this section we will look more closely at the situation in 
which the domain Q is a cube and the mesh is a uniform cubic mesh. For this case 
the dual mesh is also uniform cubic, being shifted by one half the mesh spacing in 
each coordinate direction from the primal mesh. 

For this uniform mesh case the basic error estimate can be improved. The difference 
in the uniform mesh case arises from the improved approximation theory which is 
possible. Apart from approximation theory the earlier results have obvious analogs 
which it is unnecessary to spell out in detail. The new approximation results corre- 
spond to lemma 2 and theorem 5. 


Lemma 3 Assume that u € H^(Q.) - Then we have 


€ C1) _ „(2)l 


<Ch 2 \u\ 2 #. 
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Proof. As previously, let k 3 denote the polyhedron associated with the face fij and 
the dual edge a f y By a Sobolev embedding theorem we have 

t = 1 , 2 , 

where Kj = cr' and — fij- Hence V : u — nj 1 ^ - u ^ is a bounded linear functional 

on H 2 (kj). For uniform meshes it is easy to verify that vanishes for linear 

u so it follows that 

«5 1) - U J 2) | $ C ( k )\ U \2,K } 

and the usual rescaling argument gives that the right side is 

< Ch 1/2 \u\ 2 K] . 

Then it follows that 


E l*i 1) - « 5 2) | 2 Sjh'j < C Y, h |u|L/^' 

i&j efi e 

< ch<Y. HL,- 


which gives the result. ■ 


Theorem 6 Assume that u £ then we have 

|w- u(1) |h < Ch 2 H 2 ,u- 

Proof. The proof is similar to the proof of theorem 5. ■ 


6 Numerical test. To confirm the rate of convergence given by the previous 
theorem we computed a numerical example. The domain of computation is the unit 
cube. The domain was first divided into equal small cubes with dimension hxhx h. 
Then a dual mesh was generated (dual cubes). We consider the following problem 

div u = 0 

curl u = u? 

Uy | y=0 — | z=0 — 0 

sin( l)cos( jr)cos(r) 

— 2cos(x)sin( l)cos(Y) 
cos(x)cos(j/)sin(l) 


|x=0 = 

U x | j~— i = 

u y ly=l ~ 

U z \z= l = 


9 



where u = ( u x . u y . u, ) and 


cos(x)sin(j/)sin(*) ^ 
— 2sin( x )cos( y )sin( 2 ) 
sin(.T)sin(t/)cos(a) j 


u> = 


The exact solution of this problem is 


u = 


' sin(x)cos(«/)cos(.z) 
-2cos(z)sin(y)cos(a) 
y cos(x)cos(j/)sin(a) 


Four meshes were used in the computation, h = 0.5 for the coarsest mesh and 
h = 0.0625 for the finest mesh. The results are shown in the table. 


Table 1. Errors Between u and u^. 


1 h 

0.5 

0.25 

0.125 

0.0625 

lk-u (1, Iw 

0.26d-01 

0.56d-02 

0.13d-02 

0.31d-03 


The average rate of convergence for this example is about h 21 which is slightly 
better than the rate given by the theorem. 


Conclusions. We have presented an algorithm for discretizing three dimensional 
div-curl systems in three dimensions. We proved uniqueness of the solution to 
the discretized equations and an error estimate showing first order convergence in 
general and second order for regular meshes. A basic numerical example was also 
shown. 
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